threshold and dry weight temperatures are region dependent – how well does this translate to decades in the future when regional relationships may very well change?

Initial Dataset

This code block will create threshold and dry weight factor estimates for brickman points.

Caroline’s dynamic thresholds create cfin dry weight estimates based on June and August. Therefore, we will read these months and downsample by 75% to improve runtime.

bathy <- brickman::read_brickman(scenario = "PRESENT", 
                                     vars = c("Bathy_depth", "land_mask"), 
                                     interval = "ann", 
                                     form = "stars")

sst <- brickman::read_brickman(scenario="PRESENT", 
                               vars = "SST", 
                               interval = "mon", 
                               form = "stars")

brickman <- c(bathy, slice(sst, month, 6), slice(sst, month, 8)) |>
  st_downsample(n = 1) |>
  st_as_sf(coords = c("x", "y"), crs = 4326) |>
  filter(land_mask == 1) |>
  na.omit() |>
  rename(SST_jun = SST, SST_aug = SST.1)

Mapping regions

Next, we map these brickman points to regions provided by Caroline. Points outside of the regions are omitted.

regions <- read_sf(dsn = "shapefiles/Regions_dw_vd_poly_all.shp") |>
  st_make_valid() |>
  st_transform(crs = 4326)

brickman <- mutate(brickman, region_tc = 
                     st_intersects(brickman, regions, sparse = FALSE) |>
                     apply(1, function(u) ifelse(any(u), regions$id[which(u)], NA))) |>
  na.omit()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Dry weight conversion

Next, we implement Caroline’s dry weight conversion equation to determine C. Finmarchicus dry weight multiplier at each present-day brickman point. Dry weight in is units of micrograms (\(\mu g\)).

Original conversion from Caroline:

pres |> 
  mutate(Cfin_C5C = 
           ifelse(REGION_predictions %in% c("GoM", "WSS", "ESS") & MONTH %in% c(5,6),
                  (-18.8 * T0_50)+332,
                  ifelse(REGION_predictions %in% c("neGSL","nwGSL", "sGSL", "sNL") 
                         & MONTH %in% c(7,8), (-18.8 * T0_50)+332, NA)),
         Chyp_C4C = 
           ifelse(REGION_predictions %in% c("GoM", "WSS", "ESS") & MONTH %in% c(4,5),
                  ((-4.71 * T0_50)+94)*6,
                  ifelse(REGION_predictions %in% c("neGSL","nwGSL", "sGSL", "sNL") 
                         & MONTH %in% c(5,6), ((-4.71 * T0_50)+94)*6, NA)),
         Cgla_C4C = 
           ifelse(REGION_predictions %in% c("GoM", "WSS", "ESS") & MONTH %in% c(4,5),
                  ((-4.71 * T0_50)+94)*2,
                  ifelse(REGION_predictions %in% c("neGSL", "nwGSL", "sGSL", "sNL") 
                         & MONTH %in% c(5,6),((-4.71 * T0_50)+94)*2, NA)))

Implementing for Brickman:

brickman <- brickman |>
  rowwise() |>
  mutate(idw = (-18.8 * ifelse(region_tc %in% c("GoM", "WSS", "ESS"), 
                               SST_jun, SST_aug)) + 332) |>
  ungroup()

quantile(brickman$idw, probs = c(0, .05, .50, .95, 1.00))
##        0%        5%       50%       95%      100% 
## -58.07668  10.57657  98.94075 212.36640 231.51398

Manipulating to reasonable dry weight

Expected thresholds from Caroline:

Expected dry weight values

Differences are likely due to the fact that I’m using SST instead of T0-50. Currently waiting on T0-50 values to consider, but in the meantime:

brickman <- brickman |>
  mutate(bidw = map(idw + 120, ~max(82, .x)) |> unlist())

Original Threshold

Next, we find the threshold values for the brickman points. Threshold depends on bathymetry, month, and region and is determined via the text tables provided by Caroline. We will determine threshold for the months of January, April, July, and November.

Some regions are clumped together for threshold conversions.

region_thresholds <- function(region_tc) {
  if (region_tc == "GoM") {
    "gom"
  } else if (region_tc == "WSS") {
    "wss"
  } else if (region_tc %in% c("ESS", "sNL")) {
    "ess"
  } else {
    "gsl"
  }
}

brickman <- select(brickman, -land_mask, -SST_aug, -SST_jun) |>
  mutate(Jan = 1, Apr = 4, Jul = 7, Nov = 11) |>
  rowwise() |>
  mutate(region_thresholds = region_thresholds(region_tc))

Not all months perform well in the text conversion due to lack of data, so we map to adjacent months using mapping provided by Caroline.

month_map <- matrix(c(
  11, 11, 4, 4, 6, 6, 6, 9, 9, 10, 11, 11,
  11, 11, 4, 4, 5, 6, 6, 9, 9, 10, 11, 11, 
  11, 11, 3, 5, 5, 6, 7, 8, 9, 11, 11, 11,
  1, 1, 4, 4, 4, 7, 7, 8, 8, 12, 12, 12), 
  nrow = 12,
  dimnames = list(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), c("ess", "wss", "gsl", "gom"))
)

brickman <- brickman |>
  mutate(across(Jan:Nov, ~month_map[[.x, region_thresholds]]))

brickman_mon <- brickman
thresholds <- read.table("Calanus_threshold.txt", header = TRUE) |>
  group_by(month, region_vd)

tkeys <- group_keys(thresholds)
tlist <- group_split(thresholds, .keep = FALSE) |>
  map(~as.data.table(.x) |> 
        setkey(depth))
names(tlist) <- paste(tkeys$region_vd, tkeys$month)

get_threshold <- function(mon, bdepth, region) {
  tlist[[paste(region, mon)]][data.table(bdepth), roll = "nearest"]$Min_cfin_mgm2
}
brickman <- brickman |>
  mutate(across(Jan:Nov, ~get_threshold(.x, Bathy_depth, region_thresholds)))
## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).

Resting Threshold

Using new set of thresholds from Caroline, using minimum estimates.

thresholds <- read.table("Calanus_threshold_resting.txt", header = TRUE) |>
  group_by(month, region_vd)

tkeys <- group_keys(thresholds)
tlist <- group_split(thresholds, .keep = FALSE) |>
  map(~as.data.table(.x) |> 
        setkey(depth))
names(tlist) <- paste(tkeys$region_vd, tkeys$month)

get_threshold <- function(mon, bdepth, region) {
  tlist[[paste(region, mon)]][data.table(bdepth), roll = "nearest"]$Min_cfin_mgm2
}
brickman_r <- brickman_mon |>
  mutate(across(Jan:Nov, ~get_threshold(.x, Bathy_depth, region_thresholds)))
## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).

## Warning: Removed 2 rows containing missing values (geom_bar).
## Removed 2 rows containing missing values (geom_bar).

Necessary Abundance

what is the necessary abundance based on the bandaid temperature and resting thresholds?

needed_ind_b <- brickman_r |>
  mutate(needed_ind = (Jul * 1000)/bidw)

Version expressing needed abundance as percentile relative to actual dataset.

root <- "/mnt/ecocast/projectdata/calanusclimate/src"
ae <- readr::read_csv(file.path(root, "tc_datasets/ae_regionvd_cfin.csv.gz"),
                      col_types = readr::cols())

quantile(ae$abundance, probs = c(0, .05, .5, .95, 1))
##           0%           5%          50%          95%         100% 
## 0.000000e+00 3.219821e+01 3.148560e+03 6.040407e+04 1.338762e+06

Plot considering only regions with threshold less than absolute maximum